Nonlinear Mode Decomposition: 
a noise-robust, adaptive, decomposition method 
based on the synchrosqueezed wavelet transform^ 

Dmytro latsenko, Aneta Stefanovska*, Peter V. E. McClintock 

Department of Physics, Lancaster University, Lancaster LAI 4YB, UK 



Abstract 

We present Nonlinear Mode Decomposition (NMD), a new adaptive decomposition tool for sig- 
nal analysis based on the synchrosqueezed wavelet transform (SWT). It decomposes a signal 
into its nonlinear modes, i.e. into its full oscillatory components, including all harmonics. The 
NMD procedure consists of four parts, each of which is a useful technique in its own right: 
(i) accurate adaptive curve extraction from the synchrosqueezed wavelet transform; (ii) identifi- 
cation of possible harmonics; (iii) reliable identification of the true harmonics; and (iv) recon- 
struction of the nonlinear modes from the SWT. We demonstrate the qualitative and quantitative 
superiority of NMD over the empirical mode decomposition (EMD) and ensemble empirical 
mode decomposition (EEMD) methods, and we show that NMD is noise-robust. We illustrate 
its application to a simulated signal and to a human EEG recording, obtaining excellent results 
in both cases. We point out that NMD is likely to be applicable and useful in many differ- 
ent areas of research. The necessary MATLAB codes for running NMD are freely available at 
http://www.physics.lancs.ac.uk/research/nbmphysics/diats/nmd. 
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1. Introduction 

In real life, most time series are superpositions of many different components, typically cor- 
responding to different processes occurring in the originating system. For example, blood flow 
signals, in addition to a mean flow, contain oscillatory components corresponding to heart, respi- 
ratory, myogenic, neurogenic and endothelial activities |[lf|. Because these components are mixed 
in a complicated manner, and may also be interdependent, the analysis of such signals presents a 
challenging problem. For example, we cannot extract the instantaneous phase from a multicom- 
ponent signal by application of the Hilbert transform because, if we try to do so, we will obtain 
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a mixture of the phases of all components. Thus it is usually necessary to start by decomposing 
the signal into its fundamental components, so that each of them can then be studied separately. 

How best to accomplish this decomposition is a problem of long-standing, and many differ- 
ent approaches have been proposed. The oldest is the Fourier transform (FT), which transforms 
signal from the time to the frequency domain, i.e. decomposes it into a superposition of sinu- 
soids each of constant frequency. However, even in the rare cases when the components of the 
signal are sinusoidal, they might still be nonautonomous, i.e. varying slowly in amplitude and 
frequency. Because the FT is a stationary method, it cannot encompass such variations and it 
separates them into additional sinusoidal components, complicating interpretation. One way of 
tackling the problem is to compute the FT within a moving time window, either constant or adap- 
tive, leading to a time-frequency representation of the signal. Examples include the short-time 
Fourier transform (STFT), wavelet transform (WT), and Wigner-Ville transform (WVT) 10, S]. 
Of these, we will concentrate mainly on the wavelet transform, which is the most widely used 
on account of its adaptive window size and ease of interpretation. A time-frequency representa- 
tion specifies the frequency and amplitude of a signal's components at any given time. Thus, the 
components' frequency bands can be identified, either visually or by using an automated scheme, 
after which they can be reconstructed by filtration of the WT within this frequency band. An- 
other well-known approach utilizing the discrete WT is known as subband-filtering/coding, sub- 
sequently generalized to wavelet packet decomposition (WPD) Qlltl. However, in cases where 
the frequencies of the components vary significantly in time, a single component may be decom- 
posed into two, because the determination of the frequency bands in WPD is usually based on a 
fixed sampling frequency, and thus is non-adaptive. This obviously poses a problem. Because 
the two approaches are qualitatively the same, and to avoid complications, we will consider the 
former method, based on the continuous WT. 

Other approaches to decomposition have also been proposed. Of those that work directly in 
the time domain, empirical mode decomposition (HMD) |4j] has probably been the most success- 
ful. It is applicable to a single time-series. Using it, one first identifies the local maxima and 
minima of the signal and interpolates between them to find the upper and lower envelopes. The 
mean of the envelopes is then subtracted and the resultant signal is treated as a first approxima- 
tion to the oscillatory component with a highest mean frequency. Next, the procedure is repeated 
to produce a further approximation to this component, and so on, until some stopping criterion is 
met. The final approximation will be the component with the highest frequency present. It is then 
subtracted from original signal and the whole procedure is repeated until no more components 
can be extracted. Although entirely empirical and lacking a formal mathematical basis, EMD has 
successfully been applied in many fields jsj]. Because it is based on local minima-maxima iden- 
tification, however, EMD becomes inaccurate, or fails completely, in the presence of the noise 
that is of course an unavoidable attribute of any real-life signal. To overcome this drawback, 
a modification of EMD called ensemble empirical mode decomposition (EEMD) has recently 
been proposed |0]. The idea is to add a sequence of independent realizations of white noise to 
the signal and apply EMD each time, effectively ensemble-averaging each of the components. To 
a large extent this approach solves the problem of noise, as well as the mode-mixing problem |@1. 
It should be noted that EEMD still works even when colored noise is present. However, in addi- 
tion to certain other drawbacks, EEMD still suffers from noise much more than decompositions 
based on a time-frequency representation. 

For multivariate time-series one of the most widely-used decomposition methods is principal 
components analysis (PCA) fl\\. Also known as eigenvector analysis, it explores the correlation 
structure of time-series to find those parts that are common to all signals. Given multivariate 
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time-series (e.g. brain EEG measurements from different probes), PCA assumes that all signals 
are linear mixtures of the same components, represents the time series in matrix form, and finds 
the components as orthogonal modes having maximum cross-correlation between each other. 
Thus, given signals each of length L, we represent them as an x L matrix A with a separate 
de-meaned signal in each row, calculate the correlation matrix C - A^A, and then find the signal 
components as appropriately transformed eigenvectors of C. The more strongly some component 
is represented in all signals, the larger the eigenvalue will be for the corresponding eigenvector. 
The resultant components will be maximally correlated among signals and minimally correlated 
with each other. The independent component analysis (ICA) approach is based on the 

same matrix formulation, and can be regarded as a generalization of PCA. It retains the flavor of 
PCA but, instead of finding minimally correlated components, it seeks the statistically maximally 
independent ones (e.g. those having minimal mutual information) by maximizing higher order 
cumulants (while PCA maximizes a second order one - correlation). 

Many extensions to these methods have been proposed, one of the best known being kernel 
principal component analysis (KPCA) ifioll . which extends PCA to the decomposition of non- 
linear mixtures of source signals as well. This is a useful extension, because in many cases the 
components of each signal can be nonlinearly transformed, as well as mixed in complex ways. 
Approaches are available for including nonlinearity in ICA, as well as many other modifications 
and symbioses. Nevertheless, PCA, ICA and most of their modifications work ideally only when 
the number of signals is larger than, or equal to, the number of components present. In reality, 
however, one seldom has so many signals. Often only a single time-series is available, whereas 
to apply PCA or ICA one needs to have at least two. 

In fact, PCA and ICA can still be applied for single-signal decomposition if one can create 
many signals out of the given one. This is usually accomplished by "embedding" ifTlll , i.e. by 
constructing new signals as time-shifted versions of the original one. Applied together with PCA, 
this approach is generally known as Karhunen-Loeve expansion, or singular spectrum analysis. 
Unfortunately, components obtained in this way are not unique and they may depend sensitively 
on the chosen embedding dimension and time-delay 11211 . In addition to other drawbacks il2ll . 
this makes it diflicult to obtain meaningful results in practice. We will be interested mainly in 
the decomposition of single time-series, and so PCA and ICA will be not discussed further 

For decomposition of a single signal one can use either (E)EMD, or reconstruct components 
from a time-frequency representation filtered in the chosen frequency bands. Each approach 
produces components of the form A(f) cos(0(f)), incorporating the amplitude and frequency time- 
variability that are ubiquitous in real systems. However, there are two problems that greatly 
restrict the applicability of such methods. 

The first problem is frequency-overlap. When dififerent components of the signal have fre- 
quencies that, at some time, come close to each other, the corresponding WT regions will overlap. 
This makes a full separation of components impossible and it distorts the corresponding WT co- 
efficients, making the reconstruction inaccurate. The frequency ranges of the components are 
usually relatively broad, so that such overlap occurs very readily. Although acting in the time 
domain, (E)EMD is affected by the same problem, and probably to a similar extent. One advan- 
tage of the WT is that we can improve the situation by changing the trade-off between time and 
frequency resolution, whereas in (E)EMD this is impossible. 

The second problem is that the components are usually either inherently nonlinear, or become 
so due to the process of observation, i.e. they do not represent simple sinusoids with some varying 
amplitude and frequency, but are something more complicated. In mathematical terms, let (pcU) 
be the phase of some component contained in the signal. Then the corresponding component 
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does not usually have a simple linear form like c(t) - A(t) cos((pc(t) + 4>o) (as implicitly assumed 
by (E)EMD), but a nonlinear one like c(t) = rjicpdt), t), where rj is some 27r-periodic function of 
phase: ri((p + 2ji, t) — ri{(p, f) for any t. Due to its 2ji periodicity in phase 7/((^, f) can be Taylor- 
expanded 

c(f) = r]((f>M,t) = Ai(f)cos(0c(O) + A2(f)cos(20,(f) + (^2) + AjCO cos(30,(f) + .^3) + - (1-1) 

The expansion terms proportional to sin(0c), sin(20c), sin(3<^c), ■■■ are called respectively the first 
(or main), second, third, ... harmonics of the full component. Note that, while we allow for some 
change of form with time via the explicit time-dependence of the harmonic amplitudes A,(f), 
the phase shifts (/j, should be the same if the harmonics are all related to the one component. 
However, both (E)EMD and time-frequency analysis regard harmonics as independent compo- 
nents. Thus, given a single nonlinear component in the absence of noise, EMD will usually 
decompose it into its harmonics, and filtering from a time-frequency representation will do just 
the same. Apart from the obvious logical inconsistency, such an interpretation is liable to lead 
to wrong conclusions, e.g. the apparent existence of phase-synchronized components (because 
the harmonics must obviously have perfect phase coherence). Moreover, when there are several 
nonlinear components, reconstruction by filtration of the WT in some chosen band will miss 
the harmonics of the selected component, while inadvertently picking up the higher harmonics 
of one or more of the components appearing at lower frequencies. To distinguish full nonlin- 
ear components of the form ( 11.11 ). from the only partial ones (e.g. harmonics) reconstructed by 
different methods, we will from now on refer to the former as nonlinear modes (NMs). 

The first problem can to a large extent be solved by use of the recently-developed syn- 



chrosqueezed wavelet transform (SWT) 11 1311 . It is constructed from the WT and provides a 
much more squeezed time-frequency representation, so that the different frequencies can more 
easily be separated. The main advantage of the SWT is that, due to the very good frequency 
localization, one can usually extract the desired component on its own by tracing the correspond- 
ing curve (see below). Note, that noise and frequency overlap (if present) can still affect the 
SWT amplitudes (because of being based on the WT), but extraction of only considered com- 
ponent's frequency profile, and its use with ( 18. lb . to some extent solves the problem: see Sec. 
[H In comparison with (E)EMD, SWT-based decomposition brings many advantages, e.g. having 
a rigorous mathematical basis and allowing a trade-off between time and frequency resolution 
ifisll . It is also more noise-robust. Indeed, due to the SWT's compactness, one can usually distin- 
guish oscillatory components that are close in frequency even in the presence of strong noise of 
any type. Nevertheless, it still can fail when the signal is so complicated or noisy that frequency 
curves are so corrupted or entangled to the extent that distinguishing separate one becomes im- 
possible even in the SWT. However, our experience is that, in such situations (E)EMD, as well 
as any other methods, does not yield results that are more meaningful. This seems rather hard 
to prove due to (E)EMD being purely empirical, but is based on many tests. In the same time, 
using SWT sometimes it is possible to tackle even such problem by careful adjustment of the 
time-frequency resolution tradeoff in WT to better localize components. 

Use of the SWT does not, however, solve the second problem. Thus the components extracted 
will still be harmonics of the full nonlinear mode. This problem has been partially addressed by 
Sheppard et al lfl4ll . who proposed a method for the identification of harmonics. The method 
proposed is not, however, a decomposition tool but, rather, a way of distinguishing independent 
components from the harmonics of other ones. As an additional complication, extraction of the 
amplitudes and phases for the higher harmonics can be quite tough even using the SWT, partly 
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because they are blurred by the wavelet frequency resolution, and partly because they are usually 
of lower amplitude and thus more affected by noise. 

In this paper we present a method for decomposing a signal fully into its separate nonlinear 
modes. It is capable of identifying and extracting the harmonics, including high harmonics, 
accurately, even in the presence of strong noise, leading to accurate separation and reconstruction 
of the nonlinear modes. 

2. Nomenclature 

For convenience of the reader, we now clarify the terminology that we will use throughout 
the rest of the paper, as follows - 

• Signal implies a single time-series. 

• Component is used quite generally to mean a part of a signal that can be distinguished on 
the basis of some specified criterion, e.g. having a particular form (as in Fourier analysis, 
(E)EMD and NMD), or being independent from the rest (as in ICA). The meaning will 
always be clear from the context, and often we will mean a harmonic. 

• Oscillatory component means any component of an oscillatory nature. 

• Harmonic means a term of the form A(f) cos(0(r)). 

• SWT curve, frequency curve, curve all imply a region of the synchrosqueezed wavelet 
transform that is separated from the others and contains an harmonic (as in Fig.|2]i. 

• Frequency profile or Instantaneous frequency means the extracted time-dependent fre- 
quency of a harmonic. 

• Nonlinear mode or: full nonlinear component implies a component of the form dl.lb . con- 
taining all of its (dependent) harmonics; harmonic and nonlinear mode are synonymous 
when the nonlinear mode contains only a single harmonic. 

• First harmonic is the lowest in frequency harmonic of a given nonlinear mode. This is the 
main harmonic and the phase and frequency of the full nonlinear mode are the phase and 
frequency of its first harmonic. 

• Maiyi curve is the SWT curve corresponding to the first harmonic of a nonlinear mode. 

3. Background: wavelet and synchrosqueezed wavelet transform 

Before proceeding to Nonlinear Mode Decomposition, we first review the wavelet and syn- 
chrosqueezed wavelet transform methods briefly. Given a signal s{t), its wavelet transform 
Ws{a, t) is constructed as 

WAa,t)^—\ g(u)4—)du, (3.1) 

^|a J-oo ^ a ' 

where a sets the scale and \^{u) is a chosen wavelet function. We choose it to be Morlet wavelet 

^(u) = Bf, [e'^''f'"' - e-i(2'^/«)') e-"''\ (3.2) 

where B ^ is some normalization constant. The latter was originally set to be Z?y„ = [ ^fni 1 + 
g-(2/r/o)" _ 2e^i^~''f''^~)Y^I'^, but it can in general be chosen to be almost anything; in codes we 
use Bfg - for convenience (e.g. for ( 18. II )). Nonetheless, all formulas in the paper are 

presented for the general By„ . In (13.2b . /o is the so-called central frequency, which determines the 
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relationship between the time and the scale (or frequency) resolutions of the wavelet transform: 
effectively it changes the spread of the Gaussian window. Where not otherwise specified, we use 
/o = 1. 

From ( 13.1b and ( 13.2b it is evident that to move from a time-scale ( 13.11 ) to a time-frequency 
representation of the signal one should set in ( 13.11 ) a - fo/f, where / denotes frequency. The 
wavelet transform is characterized by its resolution in time, i.e. how rapidly changes of frequency 
can be distinguished, and in frequency, i.e. how close in frequency are the components that can 
be distinguished. According to the uncertainty principle, sharp localization in time and frequency 
are mutually exclusive, obeying A,A / > c, where c is a constant dependent on the exact definition 
of At J 12. 0]. The inequality becomes an equality only for a Gaussian window 12. Si. so the 
Morlet wavelet (13.21 ) is commonly used because it maximizes the joint time-frequency resolution. 
The ratio between A, and Af can be varied by changing the central frequency /o in ( 13.21 ). which 
allows for a trade-off between time and frequency resolution. The wavelet transform does not 
lead to loss of information, so that the original signal can be reconstructed from its WT as 

/ /U — 1\ da 

dt 4 )w,{a,t)—, (3.3) 

oo J{) ^ a ' 

where the constant in front of the integral is determined as D^j, - with ip{^) denoting 

the Fourier transform of the wavelet ( 13.2b as 

Xoo 
i/r(M)e-2™'ft/M = V2^B/„ (e-2^'(A-f)' - g-Z^r'a^f')) . (3.4) 
00 

Note that the second term in ( 13.2b . which is usually small and is often ignored, nonetheless 
assures the admissibility condition, i.e. makes the Fourier transform of a Morlet wavelet zero for 
^ — > 0, which is very important, e.g. in the calculation of in (13.7b (see below). 

In real cases, both time and scale (frequency) take discrete values. Time tk - kAt - k/fs 
where is the sampling frequency of the signal. Scale values are usually chosen to be equi-log- 
spaced, and here we use the dyadic convention 

fl; = cr'^'fli s 2~fli, (3.5) 

where the so-called "number of voices" «,, - 1/ log2 cr determines the number of scales in the 
dyadic interval, and ai = 2//, is the smallest possible scale (corresponding to the Nyquist fre- 
quency, see below). All the above formulas are easily generalized to the discrete case. 

Of course, real signals are not infinite in time, so we must make a convolution between the 
signal and some part of the in-principle time-infinite wavelet function in (13.1b . This yields a 
cone of influence, i.e. a range of scales and frequencies for which the WT is determined with 
some predefined accuracy. Consider a = 1 . Then, assuming the recording to be long enough, for 
the starting time f = 0, the wavelet transform is constructed by convolution of only half of the 
wavelet function. The possible error in such a WT coefficient will be 100% because, if we look 
at the modulus of the wavelet function, the unused part (Gaussian tail for r < 0) will be equal 
to the used part. Thus, depending on the unknown behavior of the signal for f < 0, the wavelet 
amplitude obtained 1^^(1,0)1 may be twice the value found, or may drop to zero. There will 
also be a large uncertainty in the determination of its complex phase. Moreover, the mean value 
of a half-wavelet is not zero, introducing additional complications. Under the same conditions 
but for t - 2, |Wj(l,2)| is determined with » 2% precision (the ratio of the unused and used 
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parts of the wavelet modulus). It is worth noting that in the case where the signal has different 
behavior for f < 0, e.g. higher amplitudes of components, the inaccuracy might be larger The 
same considerations apply at the end of the signal. Often one "pads the data" at both ends of 
the signal to make it longer (usually by a factor of 2 to use fast Fourier transform algorithm), but 
such choices are artificial, and there is no universal padding that improves accuracy in all cases. 
We use here a precision of a; 0.15%, which corresponds to the wavelet function having at 

least three physical (overlapping with the signal) oscillations for u - t > and for m - f < 0, and 
we regard all other wavelet coefficients as unphysical. In addition, the frequency range of the 
wavelet transform is restricted from above by the Nyquist frequency f„,m- - f jl and from below 
by the lowest frequency satisfying the accuracy condition, i.e. where the wavelet scale becomes 
so big that its central part akeady contains the whole signal; for the chosen accuracy and /o = 1, 
it is roughly /j„,„ - 6/T, where T is the length in time of the overall signal. 

However, the wavelet transform frequency localization at /o = 1 is often not enough to dis- 
tinguish, even visually, the frequencies of different oscillatory components (see Fig. \lla)). We 
can of course increase /o, but then we will lose time resolution, which might lead to the sit- 
uation when oscillations in the frequency of some harmonic are in themselves regarded as an 
independent harmonics, as in the case of a Fourier transform. However, the recently introduced 



synchrosqueezed wavelet transform (SWT) 111311 is capable of providing a time-frequency rep 



resentation with much better frequency resolution and the same time resolution (see Fig. [Ub)). 
It is constructed in the following way. First we identify the frequencies f(a, t) with which the 
phase of a wavelet coefficient is growing for each scale and time: 

1 d 

f{a,t)^—-&rg{WAa,t)), (3.6) 
2n ot 

where arg(-) denotes the (unwrapped) phase of a complex coefficient (i.e. z - x + iy = 

and the multiplier 1 jln is needed to convert from circular frequency to normal. Note that in the 

original paper |13| (as well as in others like ifisil) the authors determine frequency as f{a, t) — 



-i(Ws(a, t)y^ -^^Wsia, f) (or just the imaginary part of the corresponding expression). However, 
in practice there are problems with the stability of such frequency estimates due to there being 
places with almost zero |Wj(a,f)l, where the corresponding estimate becomes singular. So one 



needs to impose some threshold on |Wj(a, t)\ 11511 . Calculation of frequency by ( 13.6b avoids such 



problems, and is more intuitive and understandable. 

When the frequencies f{a, t) have been determined, we choose frequencies /j and form 
frequency bins around them as [/)",/)^] with bin borders ff" - - '^^^ "' fi- The syn- 

chrosqueezed wavelet transform Tsifi, f) can then be calculated as 

TAfiJ)^C-^' Ws{aj,t)afl\Aa)j, (3.7) 

j:f7<]\a, .!)<]■* 

where Acij are the distances between two adjacent scales aj, and we sum over all scale indices 
for which the calculated frequency ( 13.6b lies in the /th frequency bin at time t. The constant 
in front of (13.7b is needed to give the modulus of the SWT a clear meaning as the amplitude 

of the corresponding component. It is defined d& = I; iA(^)^ with i^(^) being the Fourier 
transform of the wavelet function ( 13.4b . Since we use a logarithmic frequency scale in the wavelet 
transform, we do the same in the SWT, taking the centers of frequency bins as being the same 
as the corresponding wavelet frequencies fj - fo/a,', it seems appropriate, and gives better signal 
reconstruction than the use of different frequency scalings in the WT and the SWT based on it. 
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Figure 1: Comparisons of the of wavelet and synchrosqueezed wavelet transforms of the signal s{t) = cos(2;r? + 
0.5 sin(2n-f/10)) + 0.75 cos(3;r?) + 0.5cos(5n-f + sin(;rr/2)), corrupted by white noise with a 25% NSR (i.e. the stan- 
dard deviation of the noise is 25% of the signal's standard deviation). The signal was sampled at 50 Hz for 100 s, we 
use riy = 64, and the time-frequency region falling below the chosen accuracy is shaded uniform gray in each case. In 
the SWT three distinct frequency curves are resolved much more clearly than in the WT. Note that, although the middle 
component has a constant frequency of 1 .5 Hz, its frequency profile in the SWT exhibits a small variation due to noise 
and to overlap with nearby frequency components in the WT. 



Taking (13.5b into account, in ( 13.7b we use a'^^^'(Aa)j - iSi^aj'^^, as suggested in jlsd, which 
seems to be slightly more accurate. Using some possibilities in MatLab, such as sparse matrix 
representation, we were able to further increase the computational speed of the fast algorithm 



suggested in 111511 . making the processing time for synchrosqueezing comparable with that of 
wavelet transformation. 

Because we disregard the wavelet coefficients outside the cone of influence as being un- 
physical, we only use well-defined WT coefficients to calculate the SWT. This avoids spoiling 
the SWT coefficients of higher frequencies by inaccurately determined phases, and hence fre- 
quencies of the lower ones, whose accuracy is lower. Additionally, since the main contribution 
to Ts{fi,t) usually comes from WT coefficients of around the same frequency Wsifo/fiJ), one 
should have at least the same (or a narrower) cone of influence in the SWT as in the WT. 

The original signal can be reconstructed from the SWT simply as 1IT3I1 

s(t) = Re 2 TAfi, = 12 TAfi, t)\ cos (arg ^ TAfi, t)). (3.8) 

Note that, in the original work, the multiplier C^^' appears only in the expression for the recon- 
struction ( 13.81 ), whereas here we have moved it into the definition of the SWT ( 13.71 ). It was done 
mainly for the sake of convenience since, in addition to giving the SWT a clearer meaning, one 
does not need to know the parameters of the wavelet used (e.g. /o in ( 13.2b ) for the reconstruction 
(ESI). 

Because the SWT is a very squeezed time-frequency representation, it provides the possibility 
of extracting by itself some chosen harmonic by selection of the corresponding frequency curve 
in the SWT. Having found the frequency interval in which the harmonic resides in the SWT 
at each time (see e.g. Fig. 2 below), it can be reconstructed by summing in ( 13.8b only over 
the determined frequency indices at each time. This is equivalent to reconstruction from the 
SWT, filtered within the time-dependent frequency interval, retaining only the support of the 
corresponding curve and setting all other SWT coefficients to zero. Obviously, the reconstructed 
harmonic will be of the form A{t) cos((^(r)), like components reconstructed by (E)EMD. Thus, 



in addition to improving the time-frequency analysis, the SWT provides an effective and noise- 
robust (unlike (E)EMD) decomposition method. At the same time, curve extraction of this kind 
is frequently impossible with the WT, as is clearly seen in Fig.[T] 

4. Basic Idea of Nonlinear Mode Decomposition 

Having reviewed the basics of the WT and SWT methods, we can now proceed step-by- 
step to the construction of NMD. SWT gives us the possibility of accurately extracting single 
oscillatory components that have varying frequency and amplitude. Usually, however, these 
components are only the lower harmonics of some nonlinear mode (11.11 ). as discussed above. So 
it is highly desirable to extract a complete component and not just one of its harmonics. The 
extraction and filtering out of such nonlinear modes can assures that further extracted curves are 
not harmonics of previous ones. The procedure for full nonlinear mode extraction entails four 
distinct sub-procedures: 

(a) To extract the main curve (first harmonic), a method for accurate, adaptive, automatic, 
extraction of the curves from the SWT is needed. 

(b) Having extracted the main curve, we need a method for extracting the candidates for its 
harmonics. 

(c) We then require a test to determine whether any particular candidate is in fact a true har- 
monic of the main curve. 

(d) Finally, we need a method for accurate reconstruction of the full nonlinear mode from 
SWT main curve and its harmonics. 

These individual sub-procedures are explained in detail below, in Secs.|5]-[8] 

5. Method for curve extraction: adaptive self-learning 

The central foundation of NMD is curve extraction, i.e. finding the region in the SWT that 
contains the required component (of the form A(f) cos 0(f)). Consider the SWT TsifiJk) of a 
given signal s{t). It will have some components that are more pronounced, and some less, i.e. 
having larger or smaller amplitudes \Ts{fhtk)\- To extract one of them, we need first to find 
the frequency region that it occupies in the SWT at each time, i.e. the corresponding frequency 
curve. In general, the full SWT will have a lot of peaks (maxima in the absolute value of the 
complex Tsifi, ft)) at any time tk- Different components will be concentrated in regions related to 
the supports of these peaks. More precisely, we will regard the support of the peak in \Tsifi, tk)\ 
at a time as a set of frequency bins specified by indices {/)supp - {i\,h + l,--,'2), or the 
corresponding frequency region {/Isupp = i/i, -,/;,), for which \T,(.fi e {/Isupp, tk)\ > and 

\T s{fit-\,tk)\ - \Ts{fi^+\,tk)\ - 0, i.e. the region in the SWT with nonzero amplitude that drops 
to zero at ends. The index corresponding to the maximum \Ts{fi, fjt)| will be called the index of 
the peak, ipeak- On account of the omnipresent noise, the modulus of the SWT may not drop 
quite to zero, however, and so for the numerical implementation we use a weaker condition in 
which the support is defined as the smallest region around the peak in question, at ends of which 
the absolute value of the SWT falls below some threshold (usually 1 % of the peak value) and is 
consistently non-increasing. 
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The SWT on this support rs({!)supp, 4) contains all information about the corresponding com- 
ponent at time tk, so its value at can be reconstructed by use of ( 13.8b . So, to extract the full 
component we need to find the supports of the corresponding curve at each time tk : {i)supp. This 
is what we mean by curve extraction. Without loss of generality we can narrow our problem 
to the extraction of the main curve only, i.e. of the most prominent one, because it can then be 
filtered out, after which the next most dominant curve can be extracted and filtered out, and so 
on. So one might think in terms of selecting the support corresponding to the largest peak on 
each iteration. In reality, the different components will vary in strength, overall support, and 
other parameters, so that this approach could lead to a situation in which some extracted curves 
correspond to one component, and others to others. Hence we need some scheme for choosing 
supports which correspond only to the component needed. 

It is usually easy to select the dominant SWT curve visually but, here, we need an automated 
procedure for most of the tasks, and in particular for building a general decomposition method 
based on synchrosqueezing. However, this is less easy than it might seem, given the impres- 
sive processing power of the eye and brain. A few methods have already been proposed for 
automatic curve extraction lfl5ill6ll but they suffer from a number of drawbacks. In particular, 
they are all parametric, i.e. they include arbitrary user-defined parameters that usually depend, 
not only on the signal properties, but also on numerical parameters such as the sampling fre- 
quency. Because there are many possible curve shapes and properties, such methods are highly 
nonuniversal, requiring parameter redefinition for each type of signal and digitizing. We need an 
adaptive method: it should update its parameters during curve extraction in the light of new in- 
formation, i.e. should have an integrated "self-learning scheme", basing further curve choice on 
the information derived from the already extracted curve part. The chosen rule should converge 
to the actual curve properties in the limit of infinite time. We may expect such a method to be 
more accurate and more universal than existing approaches. 

Without loss of generality, let us reduce the problem of extracting the consequence of the 
curve supports to the consequence of the peak indices, since after that we can readily find 
their supports. Suppose we have found the peak indices for some previous times {/peaklA: = 
'peak(f*),^ = 1, - ^K, and now need to find the next one /peak(rA'+i) on the basis of the already- 
extracted information. The most appropriate and general scheme for doing this is to select the 
peak maximizing some functional conditioned on previous values 

'peak(f^:+i) = arg max{A[/, Tsifi, tg+i)] ■ F[i, Tsifi, tg+i), \i ], (5.1) 

where A[i, T^ifi, Ik+i)] determines the "power" of the peak, while F[i, Tsifi, Tat+i), {/peaklA-] is a 
functional which reflects the probability that ; is the next peak index related to the curve formed 
thus far from all the previous ones {ipeaklA". If at some step we find by ( 15.1b that /peak is not in 
fact representing the current peak (i.e. \T ,{fi^^,^^+u tK+\)\ > \Ts{fi,,.^,tK+i)\ or \T,{fi^^^utK+i)\ > 
\Ts{fij^, tK+\)\), we should follow to the nearest actual peak and assign its index to the current 
ipeak- This is the general form of the adaptive curve extraction method. 

Note that, in general, the curve's properties are defined, not only by the peaks, but by their full 
supports. Instead of peak indices one can extract whole peak supports based on the previous ones 
('peak(0 — > {suppKO in ( 15.1b ). Howcvcr, extracting the peaks places more weight on the com- 
ponents with better frequency localization, as well as substantially increasing the computational 
speed. 

One possible way of choosing F[i,Ts(fi, tx+i), {ipeak}K] is to select A^^ relevant parameters 
Pr,r = 1, ■.,Np, determined from the frequency index for the given SWT and possibly from the 
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indices of the previous peaks {/peaklA', and to rewrite the probabihty functional as 



(5.2) 



where the fl (p,-) are some chosen functions, dependant on the whole previous history {/peaklA: 
and which suppress deviations of the next parameters from the set of values taken by previous 
parameters. 

An important issue in the scheme is that of the boundary condition, i.e. we need to retain 
some freedom of choice in ( 15.1b . ( 15.21 ) when we are just starting to follow the curve. As new 
information becomes available, the functionals / ) are each updated, restricting our further 
choice according to the curve's properties. The choice of the starting form and updating 

scheme fr''\pr) — > fl^^^\pr) are clearly of crucial importance. 

In summary, we need to make four choices: 

1. Choose a power functional A [/, Tsif,, tk)] in ( 15. 11 1. 

2. Choose the parameters pr for curve characterization. 

3. Choose the respective probabihty functions fl^\pr) in ( 15.2b . 

4. Choose an initial form and updating scheme for fl^\...). 
We now consider each of these choices in more detail. 

5.1. Power functional 

We choose the power functional as the amplitude of the peak 



Other choices are also possible, however, like \Tsifi, tkW , \og{\Ts{fi, fjt)!), but ( 15.3b seems to be 
the most appropriate. 

5.2. Relevant parameters 

We choose two parameters {Np - 2) for curve characterization in ( 15.2b . namely the present 
frequency, and the frequency difference between the present and previous peaks (proportional to 
the time derivative of frequency at that point): 



The frequency value /(f) defines the region in which the frequency of the component resides, 
suppressing transitions to other curves, while the frequency difference A/(f) defines the range of 
possible frequency jumps making extracted curve continuous in frequency. One might choose 
more parameters, like the higher derivatives of frequency (e.g. the second one to make the fre- 
quency derivative continuous as well), or the peak amplitude, or its derivatives. However, the 
choice ( 15.4b seems to be sufficient for most cases. 



A[i,T,{fi,t,)]^\T,{fi,t,)\. 



(5.3) 



P\[i, Tsifi, tK+l), {/peak)/f] = /(f^+l) = fi, 

PlU, Tsifi, tK+l), {'peakl/f] = ^f(tK+l) ^ fi ' fi,,,^{r^)- 



(5.4) 
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5.3. Probability functional 

We choose the functionals /J 2 for / and Af as probabilities of the corresponding values 
already calculated in the process of curve extraction. Thus, the probability functional becomes 

F[i, TAfi, tK^l), {/peaklif ] = Pf(fi)PAf(fi " /(Ik)), (5.5) 

where Pf and Pa/ are the respective direct probabilities (% of a chosen value among those 
already encountered). In such a manner we assure that the peak frequency and corresponding 
frequency difference are characteristic of the current component and related only to it. 

Other choices are also possible, like a Gaussian distribution for Pf and an exponential for Pa/ 
based on averages and standard deviations of frequency and frequency difference. The latter can 
either be calculated directly at each step or, much faster, roughly updated assuming some form 
of the distributions (e.g. bimodal) at all the previous steps. However, the above choice of direct 
probabilities seems to be most appropriate and intuitively clear Note that, in fact, one might as 
well use the direct joint probability P(/], Afj) that the frequency and frequency difference take 
particular values at the same time (while in ( 15.51 ) we regard Pf and Pa/ as being independent) 
and, moreover, the updating scheme (see below) will be almost the same. Nonetheless, the choice 
(15.5b seems to be sufficient and can provide a Uttle more flexibility by allowing a choice between 
two independent initial distributions. 

5.4. Initial distribution and updating scheme 

Updating the probabilities P/a/ is straightforward. The frequency range is already divided 
into Nf bins (centered at fi), and so the probability P / can be represented as Nf numbers each 
corresponding to the probability that the peak frequency is in the ith bin: Pfifd - Pf{i)- As 
to the Pa/, we divide the range of possible frequency differences into 2Nf - 1 bins centered at 
±(/;-/i) (preserving a logarithmic scale of frequencies). Taking into account that / = 2''"''^"'/i, 
we wiU have bins for frequency difference centered at [-0Nf,Nf-i,...\)in,. _ \)f^^Q^ 0\X...Nf)in,. _ 
l)/i]. Given the frequency difference A/, we can find the corresponding bin as y'(A/) - Nf + 
sign(A/)floor(«,, log(l + |A/|//i)), where floor(-) denotes rounding to the lower integer 

Each time we find the peak frequency and frequency difference according to ( 15.1b . ( 15.31 ). 
(15.51 ). we calculate the corresponding frequency and frequency difference bins // a/ and add unity 
to the probabilities P/(//) and Pa/('a/) respectively, to be reused in subsequent steps. Because 
we are interested only in the relative probability, normalization is not needed. In this way, the 
probabilities soon converge to the actual distribution, assuring complete adaptation. Moreover, 
such an updating scheme needs only a few calculations at each step, making it relatively fast. 

The only remaining choice relates to the initial probability functionals P/,a/- Clearly, it is 
most appropriate to choose them to be flat at the beginning: P/(0 - qf, Pa/(0 - qAf for V/, 
unless of course we have or choose to assume some preliminary information to the contrary. The 
larger the zero level ^/a/. the longer it takes to give more weight to some bins relative to others, 
but the adaptation to the change of curve properties is improved, and vice versa. For the scheme 
to be independent of the sampling frequency, the starting values should depend on signal length 
(in samples) or on sampling frequency. Thus, for example, if qf = L/1000 where L = fsT is 
the signal length in sampling points then, after extracting the curve during first 1/1000 of signal 
duration, some frequency can at most become twice as probable. The starting level should also 
depend on the number of voices n,,, because the finer the frequency scale becomes, the less will 
be added to each bin. Quite a good choice, capable of handling most of the cases encountered 
in practice, is ^/ = and ^a/ = j^- The parameters ^/,a/ represent the only user-defined 
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parameters in the whole scheme, and they have a very clear meaning. Namely they determine the 
localization in frequency and the curve smoothness. To make the curve narrower in frequency, 
one should decrease q f, while to make it smoother j should be made smaller. Note that, in 
some cases, these might be opposite to each other. 

Other choices of updating scheme and starting distributions are also possible. Thus, each 
time we find a peak index we might add to the corresponding probability, not just unity, but a 
value based on the current probability, or a factor such as P/(/peak) - Pfiipeak) + C/Pf(ipsak) 
which would make the scheme more adaptive to changes in curve properties albeit less localized. 
The choice of initial distribution in such a case should obviously also be modified. Nonetheless, 
the proposed direct updating scheme seems to be the most appropriate. 

However, the choice of initial distribution can be modified if some prior information about 
the system is available. Thus, to prevent some frequency bin from being selected one can set the 
initial probability of the corresponding frequency to zero. At the same time, if one finds a peak 
in the time-averaged wavelet power, then to extract the corresponding component it is enough to 
center the initial frequency probability distribution around the corresponding frequency, e.g. set 
it to be Gaussian with its center at the peak frequency and a deviation based on the peak width. 
Nonetheless, a flat initial probability distribution is usually sufficient and we use that unless there 
is reason to do otherwise. 

5.5. Implementation 

In summary, we extract the peak indices by using ( 15.11 ) with ( 15.31 ), ( 15.4b and ( I5.5l l, and im- 
plement the updating scheme described above. However, there is also the question of where 
to start, since the properties of the signal might vary in time, so that the curve appearing to 
be most prominent at one time might almost vanish at another time. We therefore find those 
time and frequency indices k,„ax, imax for which the amplitude of the sum of three (not one, to 
avoid selecting some spurious feature due to noise) consecutive SWT coefficients are maximal: 
\Ts{fi,„„, fk,„„-\) + Ts{,fi„,„, tk„„,) + Ts(,fi„,„, ?*,„„,+! )l - max. Starting from here (?*:,„„, ), we follow the 
curve forward and backward in time. We stop when we reach the limit of current accuracy (the 
gray-shaded region in Fig.[T]). 

Often, that is enough. For complicated cases, however, where the curve may exhibit some 
discontinuities and other strange behavior, curve extraction by this means can be less accurate. 
To improve accuracy, after the first run one can iterate the procedure, following the curve back- 
wards and forwards, using as the initial distributions Pf,hf those obtained from the previous 
curve followed in that direction (possibly multiplied by some propagation constant) until they 
converge, i.e. the next curve is the same as the previous one up to some accuracy. On each it- 
eration, the probability functions become more appropriate, thus assuring convergence. After 
this one might find places where for subtle reasons the curves followed forwards and back- 
wards do not agree, i.e. where there exist a few different ways of tracing a curve. Such prob- 
lems can be "fixed" by selection of the path which maximizes, e.g. the sum or the product of 
A[T,(f(tu\ tk)]Pf(f{tk))PAf{f{tk) - f{tk-i)) over the corresponding fc's. 

When the dominant curve does not exist for the whole time, or where it almost vanishes at the 
beginning or end of the SWT, the extracted curve might undergo jumps to other curves. To avoid 
this happening, one can cut out those parts corresponding to the other curves. It is sufficient 
to select the biggest frequency jumps and to study the curve properties just after and before 
them. If the curve properties are completely different (e.g. ffieir 95% ranges of frequency do not 
intersect), the corresponding part is regarded as the other curve and is excised. Nevertheless, 
usually it is enough to make only one run without either of the iterative procedures, fixing and 
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SWT of ECG SWT of Respiration 




Time (s) Time (s) 



Figure 2: Examples of curves extracted from the SWTs of ECG and respiratory signals. The red lines show the boundaries 
of the extracted curve supports. The signals were each of 30 min. duration and were resampled to 80 Hz; we used = 64. 



curve cutting. Note that the whole procedure of curve extraction is closely similar to what one 
does visually when trying to find the curve by eye: one notes the most pronounced part of the 
SWT, follows along it to learn about the curve behavior, cuts out those parts that do not fit and, 
finally, concentrates attention on ambiguous regions where subtle choices are needed to decide 
which path to follow. 

After selection of the peak indices, it is straightforward to select the corresponding peak 
supports and extract a curve. By virtue of the method used for support extraction, the curve 
obtained is maximally squeezed and cannot contain other components. This is a great advantage 
over existing methods which use a constant window around the peak as a support, and where 
other curve parts entering this window can inadvertently be picked up. After extraction, the curve 
can be subtracted from the SWT and the whole procedure repeated to extract the next dominant 
curve. Note, that such a procedure is not too computationally expensive, the computational cost 
being comparable to that of the SWT transform itself. This adaptive scheme is able to extract 
the dominant frequency curves accurately from such different signals as ECG and respiration 
(measured by a belt around the thorax llTIl ') in all cases and without any changes, while existing 
methods need different and carefully-adjusted parameter sets separately for ECG and respiratory 
signals. In addition, the curves extracted by the adaptive method presented here are also much 
more accurate. An example is shown in Fig.|2] 

Having extracted a curve (i.e. the consequence of the peak supports {/)supp(fA)), one can cal- 
culate its instantaneous frequency /(f/t) amplitude A(fjt) and phase 0(fi) as 

j,^^^^ _ Zie{iU,„(t,)fi\Ts(fi,tk)\ 
Ijie{ih,„{ii,)\Ts{fi,tk)\ 

A(f..) = | ^-^^^'^'4 (5.6) 

0(fi) = arg{ TsifiJk)}- 

'e('lsijpp('i) 

Note that we calculate the instantaneous frequency from the whole curve support, and not just 
as the frequency of the peaks as is sometimes done. This is because the frequency curve is 
in many cases asymmetrically smeared around the actual frequency, so that it and the peak do 
not coincide precisely. The instantaneous frequency calculated according to (15.6b . although being 
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more empirically than mathematically derived, provides more accurate values. However, it seems 
that a more accurate and mathematically rigorous expression for the instantaneous frequency can 
be derived. 

The expressions for the amplitude and phase are exact and follow directly from the recon- 
struction formula ( 13.81 ). Thus, by curve extraction one can determine phase and frequency of even 
very complicated signals. For example, the ECG instantaneous phase and frequency - with the 
effective sampling being equal to the sampling frequency of the ECG, and not heart rate as in the 
marked events method - can easily be determined from the extracted first harmonic curve (e.g. 
as in lfl6ll : see also Sec.fTOt. This task gives rise to problems for other methods, e.g. the Hilbert 
transform, which are also more vulnerable to the effects of observational noise. Note that the 
phase extracted by use of Eq. ( 15.61 ) is in fact the true phase: for a discussion of the difference be- 
tween the protophase and the phase see llSll . This is so because invertible signal transformations 
of a single NM (when original signal can be unambiguously recovered from transformation, e.g. 
sin(0(r)) — > exp[sin(0(f))], but obviously not sin(0(f)) -> sin^(0(f)) = i[l + cos(20(r))], which 
is non-invertible) do not change the frequency profiles of harmonics. Only the amplitudes are 
changed and constant phase shifts are introduced. 



6. Extracting candidates for curve harmonics 

Once the main curve has been extracted, we need to find its harmonics. Possible h = 2, 3, ...th 
harmonics should be located at frequencies f^'^^ti^) = hf^'^\tk), where f^'^\tk) is the frequency of 
the main curve (i.e. the first harmonic), determined from ( I5.6l l. Thus, to extract /ith candidate 
harmonic one can just find the peak in \Ts{fi,tk)\ nearest to hf^^\tit) for each time and extract 
corresponding supports. 

In principle one can easily find the nearest peak, irrespective of its amplitude. However, due 
to noise and other factors the frequency profile may be distorted; in addition, small spurious 
noise peaks may occur in the vicinity, causing this approach to be inaccurate. It seems more 
appropriate, therefore, to proceed in a similar way to that used for extracting the main curve 
( 15.11 ), but using a different probability functional, now conditioned on the main curve frequency 

(h) 

at the relevant time. We choose the rule for finding the index of the peak i^^^itk) corresponding 
to a possible hth harmonic to be 

!<''\(f,) = argmax||r,(/;-,f,)|exp(-/; l-^' ~ ^^-^"'^^1 )| (6.1) 

where Risa free parameter suppressing deviations from the estimated frequency of the harmonic. 
The functional in (16. Il l is constructed in such a manner that it goes to zero when we move to the 
region of the next possible harmonic. We choose = 10 in ( 16.1b and restrict the region of search 
for 'pgaij(fyt) to a region of width equal to the main frequency. If we select as i^^^itk) not the actual 
peak, but the index in its vicinity (which may well happen if R is too large), we follow to the peak 

(h) 

and assign its index as ipj^^itk)- Next, for each peak, we find the corresponding support and thus 
extract a curve corresponding to the possible hth harmonic. Such harmonic extraction is very 
precise and works well even when the harmonic is almost completely buried in noise, as will be 
shown below. Note that, usually, the extraction is quite independent of the value of R, provided 
that R is not set too small. 
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7. Identification of true liarmonics 

Having extracted the curve for the hth harmonic candidate, we need to determine whether 
it is a true harmonic of the main curve, a curve that passes nearby, or an artifact attributable to 
noise. To do so, we use the method of surrogate-testing, and test against the null hypothesis of 
the independence between the main curve and the possible harmonic. If the frequency of the 
main curve varies in time, its harmonic should vary correspondingly, i.e. its behavior should be 
dependent on the behavior of the main curve, and the surrogate test should reject the hypothesis 
of independence. 

Given that one can easily calculate the phase of the extracted main curve and the candidate 
harmonic by use of ( 15.61 ), the best choice is cyclic phase permutation (CPP) surrogates lfl9ll . 
which work directly with phase (in lfl9ll called "permutation surrogates", but to not confuse 
with random permutation surrogates and emphasize their phase nature we call them CPP). They 
are constructed as follows. First, the original phase is divided into separate cycles, i.e. periods 
during which the phase grows from to Itt. The surrogates are then created by independently 
shuffling these phase cycles, i.e. by joining them again in a random sequence, thereby destroying 
the intercycle causality. Such surrogates are intuitively understandable, quick to create, and 
powerful. In fact, they seem to be the only realistic choice since other surrogates would need 
recalculation of the SWT and repetition of all procedures. 

To find the probability that the extracted curve for a possible hth harmonic corresponds to 
a true harmonic, we calculate the (l,/2)-phase coherence between the phase of the main curve 
0*''(fi:) and the candidate harmonic curve (ji^'^XtkY 

Pi,.(<^<'\<^<'"') = I^V^""^"'^"^-'^""<"'l, (7.1) 

where L\ is the time length (in samples) of the main curve. The phase coherence ranges from 
(no phase coherence) to 1 (perfect phase coherence). The phase of the hth harmonic should grow 
exactly h times faster than that of the main curve, so ( 17.1b should ideally be unity. In reality, 
noise, finite sampling frequency, finite frequency resolution, and possible component overlaps in 
the WT, all introduce errors, so usually it will be less than unity. We denote the value obtained 
for the original phases as pi,/,(0"*, 0"'*) — po. Next we construct Ns (we use Ns - 1000) 
independent CPP surrogates of the main phase and find (1, /i) phase coherence between 

the surrogate phase and hth harmonic's phase ps=\,...Ns = P\,h(<P^s=i Wj' probability 
measure (although mathematically not the true probability) that the extracted hth harmonic curve 
is a true harmonic of the main one is quantified by the significance of the surrogate test, i.e. 
by the relative part of surrogate phase coherences which are lower than the real one ps < pa- 
For example, if we found Ns - 1000 surrogate phase coherences ps=i,..,iooo and 792 of them 
are smaller than original value po, then the rough probability that the extracted curve represents a 
true harmonic is 79.2%. We regard the extracted curve as a harmonic if the probability calculated 
in this way is > 99%. Note that the surrogate test does not depend on the magnitude of the 
phase coherence. Thus, there might be an independent component located at the frequency of 
a possible harmonic, so that the phase coherence would be high but, because it does not fully 
adjust its phase to that of the main one, the surrogate test will reject it (e.g. see below. Table 1). 
Thus, the possibility of picking up some other genuine component that is nearby in frequency is 
eliminated. 
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Nonetheless, there is one other important issue. Because we select harmonic curves as being 
those nearest to the assumed harmonic frequency, we might introduce dependence by virtue 
of the selection procedure itself. Thus, we might sometimes just pick up the closest noise or 
parts of other components, and obtain a significant dependence of its phase on the main phase 
(which will be correct because we have introduced it by the extraction procedure). However, 
noise will have a much lower coherence with the main curve. To avoid such an effect it is 
usually sufficient to select some threshold for the original phase coherence po, like 0.2, and to 
regard the apparent /jth harmonic as attributable to noise if it is below the threshold. To select 
an appropriate threshold it is often enough to find coherence values between the main curve and 
some "phantom" harmonics with non-integer /j (like li - 3/2, 5/2, ...), which by definition cannot 
be real. The phase coherence threshold is then selected as the maximum value among those that 
pass the surrogate test. The threshold depends both on noise intensity in the candidate harmonic's 
frequency range and on the amount of distortion of the main curve (also usually caused by noise). 
So the minimum acceptable phase coherence should be determined for each NM separately. Such 
a phantom harmonics test is very effective and solves the problem. 

8. Nonlinear Mode reconstruction 

8.1. First method 

Having selected the main curve and all its true harmonic curves, we can now sum them to 
form the SWT of selected nonlinear mode (NM) ( II. lb . i.e. an SWT from which all values not 
belonging to the supports of the extracted curves have been filtered out (set to zero). From it one 
can reconstruct the full nonlinear component directly by use of ( I3.8l l. We will refer it as the first 
method of reconstruction. 

8.2. Second method 

It seems, however, that the first method method is not very accurate. The higher harmonics 
usually have much lower amplitudes than the main one, and so are much more affected by noise. 
Moreover, they might have relatively large inaccuracies in the estimated frequencies (13.61 ) due 
to noise and the logarithmic scaling of the underlying WT (so that, for higher frequencies, the 
separation measured in bins between harmonic curves becomes smaller, and they overlap with 
each other and become corrupted). All this can greatly affect both the frequency profile and the 
amplitudes of the higher harmonics, sometimes even introducing discontinuities into their SWT 
curves. Consequently, amplitudes, frequency profiles and phases of the higher harmonics are 
seriously corrupted in the SWT, making direct reconstruction by the first method quite inaccurate. 
To avoid misunderstanding, it should be noted that the signal can still be reconstructed perfectly 
from its full SWT (since the latter does not lead to loss of information). The complications arise 
only for reconstruction from frequency curves because, due to the inaccuracies of frequency 
estimation for the SWT, mentioned above, part of the harmonic can become spread over the 
whole frequency region of the SWT in a complicated manner. 

The situation can be improved by the use of a higher /o in (13.2b . but then we sacrifice time res- 
olution. A possible solution is to base reconstruction on the properties of the main curve (which 
can be found more precisely), determining the parameters of the harmonics through them. Thus, 
the frequency and phase of the /ith harmonic /*'''(f), 0^'''(f) should equal h times the frequency 
and phase of the main curve f^^\t),(p^^\t) (plus a constant phase shift in the latter case, iph in 
(II. lb ). Deviations from this law can only be due to noise and overlap with other curves. Next, 
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we need somehow to determine accurately the amplitudes A*''^(f)- Instead of calculating them 
from the SWT ( 13.8b . one might find them much more accurately from the underlying WT. If we 
have a signal s(f) = A 5(f) cos(27r/f + ip), its amplitude might be found to good accuracy from 

its WT as As(t) - I ~^ -, \Ws(fol f, t)\, as can be derived from (13. II ). ( 13. 21 ) using the same nota- 

tion. Such a method is less vulnerable to inaccuracies in frequency determination (13.6b caused 
by noise and curve-overlap, which is the main mechanism responsible for amplitude corruption 
in the SWT. In a similar manner one can estimate the phase through the WT but, unUke the am- 
plitude, it seems to be more accurately determined from the SWT using ( 15.6b . Summarizing, one 
can greatly improve reconstruction by consideration of 



2/2/"> (0 



2 



/o 



(8.1) 



where (...) denotes time average, and A^o is the most probable value of the difference mod2„{<p^''\t)- 
h(p^^\t)) (determined with the optimal number of bins Nbms ~ e" ''''*(L - l)" "^ H and L denotes 
the length of the data in samples). Wrapping of the phase difference to 2n is needed on account of 
possible phase slips. However, simply calculating the phase shift as the mean from the wrapped 
phase difference could be inappropriate in some cases. For example, when the actual phase dif- 
ference is distributed symmetrically around zero, wrapping the negative part of it will give the 
incorrect mean of tt: A0o was introduced specifically to correct such issues and provide a more 
accurate estimate of the phase shift. Preliminary phases and frequencies (on the right hand sides) 
are determined from ( 15.6b and the full NM is reconstructed as smMitk) = 2/i A^'^Kh) cos(0''''(fi)). 
This will be referred to as the second method of reconstruction. 

8.3. Third method 

In (18.1b we allow for slow changes in the form of the NM's, i.e. we allow for the fact that the 
relationship between the harmonic amplitudes may be not constant over time. Alternatively, we 
could assume that its form does not change, i.e. snmO) = A(t) Y^h bh cos{h<p{t) + ipi,), and attribute 
all fluctuations in amplitude ratios bi, to noise (in the same manner as with the phase in ( 18.11 )). 
Then we need additionally to set 

A<'''(f) = <A<''>(f)/A("(f)>A"\f), (8.2) 

after computing ( 18.11 ). The NM is then reconstructed as in the second method. This third method 
of reconstruction is more accurate for NMs of constant form. It is also more accurate in the 
case of extremely strong noise, such that precise determination of the harmonics amplitudes and 
phases is impossible. 

Having described all parts of the NMD procedure, we now illustrate the method by consider- 
ation of examples. 



9. Simulation example 

First, we consider a simulated signal whose composition is precisely known. It consists of 
two nonlinear modes (NMs), and is contaminated by noise, as shown in Fig. [3] 
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Figure 3: (a) The first and (b) the second nonfineai' modes, (c) Their sum before and (d) after corruption by noise. The 
latter represents the exemplary signal that we wish to decompose, i.e. to restore (a) and (b). The modes and the resultant 
signals are all sampled at 50 Hz for 100 s. 
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Figure 4: Intrinsic Mode Functions (IMFs) obtained from the simulated signal in Fig. [5] using the EMD and EEMD 
procedures (with the latest approach to fix the sifting number to 10 Q]). In the EEMD 2500 realizations of noise with the 
same variance as originally added to signal (50% NSR) was used. Other sifting stopping criteria Ql or choices of noise 
variance do not lead to qualitatively different results. 
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Figure 5: Central parts of the WT and SWT transfomis of the simulated signal shown in Fig. [5] Hamionics of the first 
NM are located around 0.5, 1.5 and 3.5 Hz, and of the second NM at 1, 2 and 5 Hz (these frequencies are indicated by 
tick-marks in both panels). 



Fig Hlpresents the results of the EMD and EEMD procedures applied to this signal. As can be 
seen, each of them fails completely, producing neither the full nonlinear modes nor the separate 
harmonics. This is mainly due to the presence of noise, which is of course ubiquitous in real 
systems. In the absence of noise, EMD will separate the considered signal into its harmonics, 
and EEMD will do the same albeit less accurately. Thus, even in the unrealistic situation without 
any noise, EMD might lead to the misidentification of harmonics as independent components. 

Fig- H shows the WT and SWT of the signal to be investigated. As can be seen, only two 
curves (among six), corresponding to first harmonics of the two NMs, are clearly visible in the 
WT while, in the SWT, the harmonics at 1.5 and 3.5 Hz are also distinguishable. Nonetheless, 
even in the SWT, the harmonics centered at 2 and especially at 5 Hz are hardly visible; but, as 
we will see, they can nonetheless still be extracted by our method. 

As is akeady clear, the example selected is a tough one, chosen to demonstrate the power of 
the method. It includes the following complications: 

1. Both nonlinear modes include higher harmonics (the 5th and 7th), which are not well 
expressed in the SWT (see discussion of the second reconstruction method). 

2. Almost coinciding with the possible harmonics of the first NM (2th, 4th and 10th) are 
harmonics of the second one. Thus, the main component of the second NM has maximum 
frequency difference of only 0.0075 Hz from the possible (although in fact absent) second 
harmonic of the first NM, making it resemble a true harmonic. Consequently, the resultant 
signal before contamination with noise, resembles a single well-defined NM rather than 
the two that exist in reahty, as illustrated in Fig.|3] 

3. The noise intensity is x 0.49 and represents: a 50% NSR (noise to signal ratio, defined 
as the standard deviation of the noise divided by that of the signal) for the overall signal; 
a 61% NSR for first NM; and a 87% NSR for the second one. Moreover, it is a white 
noise, with contributions at all frequencies, highly corrupting the already corrupted high 
harmonics, many of which have amplitudes much lower than that of the noise. 

4. The third harmonic of the first NM (1.5 Hz) and the second harmonic of the second NM (2 
Hz) considerably overlap in the WT, causing distortion in the frequency profiles and some 
overlap in the SWT, as shown in Fig.|5] 
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Figure 6: The parts of the signal's SWT coiTesponding to the first and second NMs, extracted by the harmonic finding 
procedure. The original signal, and the first and second NMs contained in it are shown in Fig. [3] and the full signal's 
SWT in Fig.[5] 



All this makes the decomposition of such a signal a highly non-trivial task. Note that the fre- 
quencies of the components have only small variations, so that they can be readily distinguished, 
even in the Fourier transform. However, signals with more pronounced nonautonomous proper- 
ties are actually easier to deal with using NMD, because the differences between true and false 
harmonics then become more obvious. The choice of small frequency variation was deUberate, 
in order to complicate our task. 

Let us now apply NMD. We first extract the dominant curve, find its harmonics, and filter 
all NM from SWT; then extract the next curve; and so on. The two main curves are quite well- 
behaved, so their accurate extraction is no problem at all, even for less sophisticated methods. 
Next we need to extract the candidate harmonics, which is also straightforward. By a phantom 
harmonics test we found the phase coherence threshold to be 0.44 for the first NM and 0.07 
for the second one. Bearing this in mind, we apply a surrogate test to determine which are the 
true harmonics. Fig. |6] shows the resultant curves of the true harmonics composing each NM, 
as extracted from signal SWT (Fig. |5jright)). As can be seen, the nonlinear components have 
been completely separated. Having extracted the true harmonics, it is now straightforward to 
reconstruct the NMs using the chosen reconstruction method. 

Table 1 summarizes the whole NMD procedure. It shows the results of the surrogate test, 
as well as the mean phase shifts {(p^''Ht) - h(p^^\t)} and amplitude relations ( ^,,||'| ), determined 
using the third reconstruction method, in comparison with the true values. It is evident that we 
can easily distinguish true from false harmonics. Moreover, even in the presence of so many 
complications and the extremely strong noise, we are able to determine accurately the amplitude 
relations and phase shifts for both the first and second NMs, including even the higher and 
completely blurred harmonics. 



Fig.|2]shows the real underlying NMs and their reconstructed versions, using the first, second 
and third methods. Remarkably, even the second NM is quite well reconstructed. To the best of 
our knowledge, no earlier method has been proposed that is capable of accurate reconstruction of 
the even harmonics of such a signal. In the case considered, the third method gives the best result, 
because the NMs are of constant form. Where there is some variation of form and the noise is 
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Table 1 : Use of the surrogate test to determine which are the true harmonics of each of the two NMs of the signal in 
Fig. |3]d). The mean phase shifts and amplitude relations, determined by the third reconstruction method, are listed in 
comparison with their actual values. In each case, the columns left-to-right provide for the Ath harmonic: the suiTogate 
significance level; the phase coherence p; the actual amplitude ratio; the extracted mean amplitude ratio; the actual 
phase shift; and the extracted mean phase shift. The true harmonics are emboldened, italic indicates harmonics of the 
other component in each case, and " — "s denote hai'monics determined to be false. The significance threshold for true 
harmonics is 99%, and the minimum phase coherence has been determined to be 0.44 for the first NM and 0.07 for the 
second one. 
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Figure 7: The first and second nonlinear modes reconstructed by different methods (see section about reconstruction). 
The original components are shown in the upper plots and by white lines on a gray background. 
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not too strong, it might be better to use the second method. However, the second method treats 
noise effects as form variations and thus is less noise-robust compared with the third method. 

Direct reconstruction by the first method usually gives the worst results, and is not recom- 
mended. Note also, that the amplitude of the quite prominent 1 Hz component, reconstructed 
from the SWT curve, is often much lower due to overlap between 1 Hz and the tails of the 1 .5 Hz 
component in the WT (see Fig.|4|. Where they at least partially overlap, the estimated frequencies 
( 13.61 ) become corrupted, which results in the assignment of the corresponding term (containing 
part of the 1 Hz component's amplitude) to the other frequency in the SWT. Determination of 
the amplitudes from the WT, as in ( 18. lb . avoids this problem. 



10. Real life example - removing cardiac artifacts from a human EEG signal 

After demonstrating its advantages in the analysis of simulated data, we now apply NMD to 
real data. We will illustrate how the NMD procedure can filter out the components corresponding 
to one signal from another one, taking as an example the removal of cardiac artifacts from a 
human electroencephalogram (EEG) signal. In measuring the EEG, one inevitably encounters 
some cardiac artifacts, due both to the direct pick up of heart electrical activity and also to blood 



flow pulsations underneath the probes (the so-called ballistocardiogram (BCG) artifacts 112 111 ). 
The latter are particularly prominent when the EEG is being measured in a magnetic field (e.g. 
during MRI) |231, in which case they are usually filtered by ICA |f2l']. 



In our case II22I1 . however, the cardiac artifacts are quite small in amplitude, making them 
harder to distinguish. Nonetheless, they can be prominent in the frequency domain, affecting the 
corresponding measures. Thus in Fig. jSja) and (b), which show respectively the WT and SWT 
of the EEG signal, the cardiac harmonics are clearly seen. We will refer to the cardiac contri- 
bution to the EEG as its ECG image, due to the apparent "mirroring" of the ECG SWT into the 
EEG SWT, as seen in Fig. [HJb) and (c). Because of the spatially localised generation of elec- 
tric potential difference by the heart, and because of the different underlying microvasculature 
(affecting BCG artifacts), the ECG image in the EEG signal takes different forms for different 
probe positions. Consequently, given its assumption of linear mixing, ICA is unable to remove 
the image in our case, as was numerically verified by applying ICA algorithm @] to the EEG 
measurements from all four probes plus ECG signal. (E)EMD also fails. Moreover, since the 
EEG contains many oscillatory components, even NMD curve extraction becomes less accurate 
due to the frequent overlaps. 

Nonetheless, to help improve extraction, we can make use of the following procedure. We 
recall that the SWT frequency profile of the main curve does not change if we apply an invertible 
transformation to the NM. Thus, such a profile should not depend on the measurement apparatus 
(assuming an invertible measurement function). It means that the frequency curve correspond- 
ing to the first harmonic of the ECG in the EEG's SWT should have the same instantaneous 
frequency as the main curve in the SWT of the original ECG signal. In reality, however, the 
curve corresponding to the first harmonic of the ECG image is more severely distorted in the 
EEG (cf. Fig. [HJb) and (c)). Therefore, instead of extracting the corresponding curve from the 
EEG SWT, one can extract the main curve from the directly measured ECG and use its unspoiled 
phase and frequency profile for extraction of the harmonics of cardiac activity from the EEG (for 
h = 1,2, sequentially). After distinguishing the true harmonics, the ECG image can then be 
reconstructed using either the second or third reconstruction method. Obviously, for all this to 
be valid, the ECG should be measured simultaneously with the EEG. 
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Figure 8: (a) Wavelet transform of the EEG signal, (b) synchrosqueezed wavelet transform of the EEG signal, (c) 
Synchrosqueezed wavelet transform of the ECG, shown for comparison. The ECG harmonics in the EEG are clearly 
seen by comparison of (c) with (a) and (b). (d) ECG artifact in the EEG, reconstructed by the third method: the black line 
shows the ECG image and the gray one shows the pure ECG at the same time (scaled in amplitude in proportion of their 
first harmonics), (e) SWT of the EEG filtered from cardiac artifacts (ECG image). As can be seen, it is free from the 
harmonics of cardiac activity(cf. b)). (f) SWT of the reconstructed ECG image shown in (d) (not from the parts extracted 
from EEG SWT, as in Fig. 6, which are quahtatively the same although noisier). Thus the ECG haiTnonics visible in (b) 
have been clearly separated into the extracted ECG image, (g) EEG with (black Une) and without (white hne) cardiac 
artifacts. The time segment includes around 10 cardiac oscillations. The full signals are of duration 20 min. and are 
sampled at 80 Hz. In computing the WT and SWT we used = 64. We also used /o = 3 to improve the visibihty of 
cardiac harmonics, since at /d = 1 they are not so clearly seen. The EEG was measured under anaesthesia (23l . but the 
same artifacts arise under all conditions. 



In the third reconstruction method it may also be better to use for the "driving" amplitude 
(A"*(r) in ( 18.21 )) the original amplitude of the ECG first harmonic, because it is much less dis- 
torted. However, this can be done only if there is expected to be proportionality between the 
power of the ECG signal and its image in the EEG signal, and on the assumption of a static 
(time-independent) measurement function. Because there is no reason to think that measurement 
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of the EEG might transform heart activity as picked-up at different times in different ways, we use 
in reconstruction the main ampUtude determined from the ECG. By this we make the assump- 
tion of linear proportionaHty between the first harmonics of the ECG and its image. Of course, 
this is an approximation since for nonhnear transformations of the time-dependent ampHtude it 
will not be the case. Nor it is so for BCG artifacts. Nevertheless, since heart activity is quite 
stable, deviations from the proportionality law are expected to be smaller than distortions caused 
by noise in the EEG. The approach still works even in the challenging case of a time-varying 
measurement function like SffMit) — > (1 +0.1 sma)t)s^^^(t), but in that case one must use the 
amplitude determined from the image's main curve. 

In practice, we extract the ECG first harmonic (narrowing the region of search to 0.5-2 Hz, as 
in Fig.|2]i and find the cardiac frequency profile and phase using ( 15.6b . Then we use the frequency 
profile obtained to extract /j = 1, .., 20 possible ECG harmonics from the SWT of the EEG. By the 
extraction of phantom harmonics with li - 1 /2, 3/2, ... 39/2 we find the threshold for the phase 
coherence value to be a; 0.15. Bearing this in mind, we apply the surrogate test to distinguish 
the true harmonics. In the example under consideration we were able to distinguish 6 true ECG 
harmonics. We assumed a more-or-less constant waveform and applied the third method of 
reconstruction because there are too many overlaps and badly corrupted components (e.g. 4th, 
5th and 6th harmonics), for reconstruction by the second method to be stable (cf. Fig. |7]l. To 
avoid amplitude overestimation of the badly-blurred harmonics 4-6 (see the beginning of the next 
section) we estimate them from the SWT according to ( 15.61 ). The reconstructed component, and 
its SWT, are shown in Fig.[8jd) and (f) respectively. The SWT of the reconstructed ECG image is 
in close agreement with the SWT of the ECG itself, but the harmonics have different amplitudes 
and phase shifts due to the different kinds of possible transformations (even the main curves being 
slightly shifted in phase). Subtracting the reconstructed ECG image from the original EEG, we 
no longer see any cardiac harmonics. Nor are they present in the filtered EEG SWT (Fig.[8je)), 
nor in the WT or time-averaged wavelet power (not shown), assuring us that the corresponding 
components have successfully been filtered out. 

In the present example, the ECG image has a completely different form from the directly- 
measured ECG, suggesting a large contribution of BCG artifacts. The saw-tooth form has also 
been observed in other EEGs where cardiac artifacts are very prominent. However, cardiac 
artifacts depend on the positions of the EEG electrodes. Often they are much smaller and can 
hardly be seen in the EEG (S)WT. The latter case is not so illustrative, but NMD works even 
there. In this case the ECG image has a sharper peak, more closely resembling that in a directly 
measured ECG. 

Based on our EEG example, it would seem that NMD can be applied widely to cases where a 
signal contains components related to another signal. For example, one can extract both the heart 
and respiratory components, accurately, from blood flow oscillations. There are also many other 
possibilities. Note that, if there is no image of an activity, then all harmonics including the main 
one will be regarded as false. Thus, applying our method to simultaneously recorded ECG and 
respiration signals, we have not found any image of respiratory activity in the EEG, implying 
that the measurement process is almost unaffected by the mechanical action of breathing. 

11. Possible issues 

The method for extraction of candidate harmonics and identification of true harmonics is very 
powerful and accurate, allowing for distinguishing harmonics that are almost completely buried 
in noise. However, estimation of the amphtudes of such corrupted harmonics according to ( 18. Il l 
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is not so accurate (e.g. see Table 1, 5th harmonic of the second NM). This is not surprising since, 
if the wavelet power of the noise is comparable to that of the harmonic itself at the relevant fre- 
quency, it will contribute substantially to the obtained amplitude, leading to its overestimation. 
In this situation reconstruction from the SWT via ( 13.8b . which usually underestimates amplitude, 
might sometimes be less imprecise (and so we used it for the badly -blurred 4-6 harmonics in the 
previous example), but neither method gives accurate results in such a tough case. Nonetheless, 
such harmonics have low amplitudes and so do not make a large contribution to the NM. More- 
over, this problem in fact can be solved by the use of iterative scheme, which will be discussed 
elsewhere. It should be noted, that even for buried-in-noise harmonics estimation of the phase 
shift according to ( 18. II ) still works fine. 

The choice of central frequency /o in the (S)WT is of particular importance. Usually /o - 1 
is sufficient and we used it throughout the paper, except for the last example, where we have 
used /o = 3 for better harmonics localization and visibility in (S)WT. In almost every case, 
however, careful adjustment of /o in ( I3.2l i can improve reconstruction, sometimes tremendously. 
Appropriate selection of /o is really a separate topic and so is not discussed here. Sometimes, 
an increase of can improve reconstruction slightly, mainly by improving the extraction of the 
main curve and harmonics. 

In the case a few of the signals contain related oscillatory components (e.g. last example), the 
frequency profiles of related components, while having the same form, can sometimes be shifted 
in time relative to each other in different signals. For example, the cardiac pulse in the blood 
flow needs some propagation time to reach the measurement site. Nonetheless, unless the time 
delay greatly exceeds the characteristic time of the frequency variations, the extraction of candi- 
date harmonics remains accurate. This is because we are selecting harmonics in some vicinity 
according to ( 16. 11 1. Even if a large time-delay is expected, one can decrease R in ( 16. U for the 
first run. However, time-delay might affect the surrogate test, sometimes leading to a true har- 
monic being regarded as false. Therefore, if one expects a non-neghgible delay, then for better 
harmonic extraction and reconstruction, one should first find the time delay At (which in itself 
provides valuable information). This can be done by preliminary extraction of the correspond- 
ing curve (first harmonic candidate) in the signal under consideration, still using the accurate 
reference frequency profile; the time delay can then be found as the time of maximum cross- 
correlation between the extracted and reference frequencies. Due to the harmonic extraction 
procedure the estimated time-delay is not exact in the presence of noise, but the inaccuracies are 
expected to be small. Then driving parameters can be tuned accordingly, i.e. f^^\t) — > /^''(f- Af), 
0^''(f) -» (p^'^\t-At) and, if we use reference amplitude, "(f) -» A*"(f-Af). Using the adjusted 
driving parameters, the whole procedure can then be applied with full rigor. In the case of the 
ECG image in the EEG, we estimated the time of the response to be at an 80 Hz sampling rate, 
so cardiac activity is picked almost instantaneously. 

For simplicity, we also assumed that the first harmonic has the largest amplitude, so that the 
extracted dominant curve will correspond to it. In general, however, this need not necessarily 
be true, e.g. the ECG third and fourth harmonics as seen in the EEG signal can be stronger than 
the first (Fig. [8]). Such cases are relatively rare but, if such a possibility is thought to exist, one 
should check the 1 /2, 1 /3, ... subharmonics to confirm that the extracted curve is really the first 
harmonic. It is interesting to note that one can consider the inverse problem in the sense of 
identifying the main curve through its harmonics. For example, let the main component be very 
low in frequency, so that only a small part of it is within the WT cone of influence (Fig. [T]l- If 
we are able to determine from this small part at least one of its true harmonics together with the 
amplitude relation and phase shift, we then will be able to reconstruct the lower harmonic outside 
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the cone of influence using the amphtude and phase of the higher one. Of course, calculating the 
amplitude of the lower harmonic in this way is reasonable only for a component of constant form, 
so that the amplitude relation is constant, while the phase obtained always remains applicable. 

The last issue to be noted is the well-known boundary problem. Because determination of 
the WT is worse near the boundaries of a signal (Fig.[T]), we should either reconstruct NMs only 
within a cone of influence, or accept that they are less accurate near boundaries. It is common to 
all decomposition methods, and seems to be unavoidable, although some improvements can be 
made in estimation of the WT at the ends, e.g. by some specific padding. It is important mainly 
for the lowest frequencies, for which many of the WT coefficients can only be determined with 
reduced accuracy. 

12. Summary and conclusions 

Our new decomposition method (IS] consists of four individually novel parts: accurate adap- 
tive curve extraction; extraction of possible harmonics; identification of true harmonics; and, fi- 
nally, reconstruction of the full nonlinear mode. It should be noted that reconstruction of a single 
NM is not computationally expensive. Thus, the whole procedure including the extraction of say 
10 possible harmonics, often take less time than e.g. EEMD with an ensemble size of only 100. 
After the first nonlinear mode has been reconstructed, we filter out those parts of the SWT cor- 
responding to it, and re-apply the whole procedure to extract the next NM. The number of NMs 
needing to be extracted usually becomes obvious during processing, but one can devise auto- 
matic criteria to cease extracting further NMs at the appropriate stage. Since NMD extracts NMs 
in order, from more prominent to less prominent, a good choice is to stop when the power of the 
next main curve drops below some threshold and it does not have true harmonics. Summarizing, 
NMD algorithm is therefore as follows: 

1. Calculate the SWT of the signal under consideration (Section 2). 

2. Extract the dominant frequency curve from the SWT (Section 4). 

3. Extract the candidates for possible true harmonics (Section 5). 

4. Use the surrogate test and the phantom harmonics test to identify the true harmonics (Sec- 
tion 6). 

5. Using all of the extracted harmonics, reconstruct the NM (Section 7). 

6. Filter the extracted NM from the signal, either by subtracting the corresponding curves 
from signal's SWT, or by subtracting the fully reconstructed NM from the signal itself (in 
the latter case, then recalculate the SWT for the filtered signal). 

7. Repeat steps 2-6 for as many times as are needed to process the distinct NMs expected to 
be contained within the signal, or until some stopping criterion is met. 

Another use of NMD, which we demonstrated on the EEG signal, is to extract the full com- 
ponent corresponding to some reference signal from which its accurate frequency profile and 
reference phase can be derived. Such problems arise quite often, and in many scientific fields. 
Because the frequency profile of the NM should not change in an invertible transformation which 
might arise due to the different measurement processes or transformations in the environment, 
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the frequency of the needed component extracted from a reference signal should be the same 
as in the investigated one (if corresponding component is present there), although possibly with 
some time delay. With this in mind, one can extract the required frequency profile and phase 
from the reference signal, find the possible time shift, adjust the frequency profile accordingly, 
and then use it to extract the corresponding component even for badly noise-corrupted signals. 
Thus, one can apply NMD whenever the component's frequency profile and phase can be found. 
It does not matter whether this information comes from the signal itself or from another source. 
If there is expected to be time-independent proportionality between the reference component and 
its image, it may be advisable to use the amplitude of the reference component as the driver in 
order to obtain better reconstruction. 

Unlike previous methods, NMD reconstructs the full underlying components of the signal, 
including all its harmonics. It is also extremely noise-robust, and we were able to obtain accurate 
reconstructions even for a 50% NSR. Application of NMD to both a simulated example and 
to a real-life EEG recording were used to demonstrate its power. The main requirement is the 
possibility of finding the first harmonic's frequency and phase. We showed that this can be 
done either directly from the signal under consideration, e.g. by SWT curve extraction, or from 
another signal containing the same component. Because the SWT is highly compressed, it will 
normally have at least some visible components even in the presence of extremely strong noise, 
as we saw in our simulation example, so that curve extraction is usually feasible. Having found 
the component's frequency profile and phase, NMD is capable of identifying and extracting the 
harmonics, even where they are almost completely buried in noise. 

In conclusion, NMD provides a powerful, noise-robust and accurate decomposition method 
of a new type, enabling the recovery of full components (including all harmonics). Its area of 
applicability is extremely wide. Thus, it can be used wherever (E)EMD has been applied (e.g. 
references in Q |^). By providing for reconstruction of the full nonlinear modes, it promises 
deeper insights into the underlying phenomena. Moreover, its extreme noise-robustness and 
other possible uses give it a wider range of applicability than (E)EMD and other decomposition 
methods. In particular, it can be applied to almost all multicomponent signals of the kind that are 
ubiquitous in the hfe sciences, climate studies and astrophysics. 
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